Strong Resilience of Topological Codes to Depolarization 
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The inevitable presence of decoherence effects in systems suitable for quantum computation ne- 
cessitates effective error-correction schemes to protect information from noise. We compute the 
stability of the toric code to depolarization by mapping the quantum problem onto a classical dis- 
ordered eight-vertex Ising model. By studying the stability of the related ferromagnetic phase both 
via large-scale Monte Carlo simulations and via the duality method, we are able to demonstrate 
an increased error threshold of 18.9(3)% when noise correlations are taken into account. Remark- 
ably, this agrees within error bars with the result for a different class of codes — topological color 
codes — where the mapping yields interesting new types of interacting eight- vertex models. 



I. INTRODUCTION 

Moore's law has accurately described the speedup of 
current computer technologies for half a century, yet this 
speedup is slowly coming to an end due to transistor 
limitations. A promising alternative is given by quan- 
tum computers. However, the qubit manipulations re- 
quired for information processing and communication are 
prone to errors because qubits are more sensitive to noise 
than their classical counterparts. Consequently, protect- 
ing qubits has become an issue of paramount importance 
for the success of quantum computation. The effects 
of single-qubit operations can be decomposed into three 
processes, bit flips, phase flips, as well as a combination 
thereof, which can be represented by the three Pauli ma- 
trices cr x , <j 2 , and (7^ , respectively. This is in contrast to 
classical bits, which can only suffer from a single type of 
error, namely bit flips. 

More generally, the notion of a noisy channel is instru- 
mental in characterizing the disturbing effects on phys- 
ical qubits. Such a quantum channel can be described 
by specifying the probability (or "qubit error rate" ) p for 
each of the aforementioned noise types. For instance, if 
only cr x occurs, then we have a bit-flip channel. Here we 
are interested in channels of the form: 

V p (p) = (l-p)p+ ]T Pvi tr w ptr w , (1) 

w—x,y,z 

where the density matrix p fully describes the quantum 
state and the probability for each type of error to occur 
is p w G [0, 1] with p := p x + p y + p z . The depolarizing 
channel exhibits equal probability p w = p/3 for each er- 
ror type to arise. As such, the depolarizing channel is 
more general than the bit-flip channel, because it allows 
for the unified, correlated effect of all three basic types 
of errors. The implications of this error model for the 
performance of a quantum code remains an open prob- 



lem. In addition, the depolarizing channel plays a fun- 
damental role in quantum-information protocols where 
noise has to be taken into account, including quantum 
cryptography p], Q, quantum distillation of entangle- 
ment [3|, and even quantum teleportation Q. Exper- 
imentally, controllable depolarization has been realized 
recently in photonic quantum-information channels [E[ 
with a Ti:sapphire pulsed laser and nonlinear crystals, 
as well as 2-qubit Bell states [6]. Here we compute the 
effects of depolarization on a set of entangled qubits. 



A. Topological codes 

The goal of quantum error correction @, [1| is to protect 
quantum information from decoherence. One approach 
using topology is based on encoding (few) logical qubits 
in a particular state subspace of (typically many) physical 
qubits which is not disturbed directly by noise. Such a 
suitable subspace of states can be defined in terms of a 
set of commuting observables, called check operators, 

^ = (tV 2 -^, (2) 

each being a projective measurement with respect to the 
code subspace (i.e., the eigenvalue signals errors on par- 
ticipating qubits). Investigating all stabilizers S l allows 
one to limit the set of possible errors to those compati- 
ble with the measured error syndrome. Our best strat- 
egy then is to classify the remaining, nondistinguishable 
errors according to their effect on the encoded logical 
information and undo the effects of the most probable 
equivalence class E. 

A hallmark of topological quantum error-correction 
codes i s the geometrical locality of these check 

operators: Physical qubits are placed on a lattice and 
check operators depend only on a few neighboring sites. 
The logical information, which is encoded globally in a 
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subspace of all physical qubits, is preserved as long as 
we can successfully detect and correct local errors. If 
errors on the physical qubits occur with a probability 
p, the error threshold p c — a key figure of merit of any 
quantum code — is defined as the maximum error proba- 
bility p, such that error classification is achievable. For 
error rates larger than p c the error syndrome's complex- 
ity inhibits unambiguous error recovery. It is therefore 
of current interest to find codes where p c is large. 



B. Error threshold as a phase transition 

The process of error correction resembles a phase 
transition and, indeed, it is possible to connect error 
correction directly to an order-disorder phase transi- 
tion in a suitable classical statistical-mechanical model 
[l2[ HEl fl6| . One can derive a Hamiltonian He of in- 
teracting Ising spins s^, labeled by a Pauli error E that 
controls the sign of the couplings, such that the proba- 
bility of each equivalence class E is proportional to the 
partition function 

p(E) oc Z E (P) := ^ HeM • (3) 

{Si} 

Equation (j3]) has to be interpreted as describing a ran- 
dom statistical model with quenched couplings and two 
parameters: the error probability p governing the fraction 
of negative interaction constants J a G {±1}, and the in- 
verse temperature f3 = 1/T. For low enough T and p 
the system orders into a ferromagnetic state (see Fig. [I]). 
Along the Nishimori line [17] where Eq. (|3)) holds, the 
ordered [disordered] phase corresponds to the topologi- 
cal code being effective [ineffective]. The intersection of 
the Nishimori line and the phase boundary identifies the 
error threshold p c . 

The first topological codes studied were toric codes ^ , 
still under intense investigation and scrutiny mainly due 
to their simplicity and elegance. To determine their error 
threshold, we show that toric codes under the depolariz- 
ing channel connect to the celebrated eight- vertex model 
(see Fig. [2]) introduced by Sutherland [l8|, as well as 
Fan and Wu [l9[, and whose general solution by Baxter 
[20l-[22j stands up as the culmination of a series of break- 
throughs in the theory of phase transitions and critical 
phenomena. 

The aforementioned mapping onto a statistical- 
mechanical model to compute the error tolerance of quan- 
tum codes was first applied to toric codes with bit-flip 
errors [l5|, connecting them to the random-bond Ising 
model. In general, for individual bit flips the error thresh- 
old is p c « 10.9%, and the same is true for phase flips 
alone. Therefore, under depolarizing noise and correct- 
ing separately bit flips and phase flips, the threshold is 
Pc — (3/2)p c « 16.4%. However, this result neglects cor- 
relations of bit flips and phase flips. We estimate the 
threshold under depolarizing noise for ideal error correc- 
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Qubit Error Rate p 

FIG. 1: (Color online) Phase boundary estimated from Monte 
Carlo simulations for the estimation of the error threshold of 
the toric code, as well as two realizations of color codes (see 
text). The error threshold p c corresponds to the point where 
the Nishimori lines intersects the phase boundary. Remark- 
ably, the phase boundaries for all three codes agree within er- 
ror bars. The stable ordered phase corresponds to the regime 
where quantum error correction is feasible. 




FIG. 2: (Color online) When computing the stability of the 
toric code to depolarization, the problem maps onto a clas- 
sical statistical Ising model on two stacked square lattices. 
In addition to the standard two-body interactions for both 
top (a) and bottom (b) layers, the resulting Hamiltonian also 
includes four body terms (c) that introduce correlations be- 
tween the layers. 



tion, such that, in particular, correlations are taken into 
account. We find p c = 18.9(3)%. Remarkably, the error 
threshold increases significantly by taking correlation ef- 
fects into account. They should thus not be neglected by 
recovery algorithms. A recent advance in this regard is 
the renormalization approach of Duclos et al. [23| where 
p c ~ 16.4% was confirmed, still leaving room for further 
improvement [24]. Note also that p c is very close to the 
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hashing bound p ~ 0.1893 [25], which is also the case for 
uncorrected bit-flip and phase- flip noise [HI, [26| . 



II. TOPOLOGICAL STABILIZER CODES 
A. Error correction in stabilizer codes 

Both toric codes [9] and color codes [10] are topological 
stabilizer codes. A stabilizer code is described by a set of 
check operators Si in the Pauli group. That is, they are 
tensor products of Pauli operators cr x , cr y , and <j z . These 
check operators Si are a commuting set of observables 
with eigenvalue ±1 that generates an Abelian group S := 
(Si) that does not contain —1, called the stabilizer group. 
Encoded states \ijj) are those for which all check operators 
satisfy Si\ip) = + If errors affect the state, typically 
they will change the value of the check operators leaving a 
trace that can be used to recover the original state. Note 
that some errors are undetectable because they commute 
with all check operators and thus leave no trace. 

We are interested in noisy channels of the form 

Po^Pi = ^p(E)E Po E^ (4) 

E 

where the sum is over all Pauli group elements E, and 
p(E) denotes the probability for E to occur. Several dif- 
ferent Pauli errors E have the same effect on the encoded 
state. Therefore, it is convenient to place them in equiv- 
alence classes E, such that E is equivalent to E' when 
EpE^ E'pE'^ on an encoded state p or, equivalently, 
when EE' is proportional to a product of check opera- 
tors. Therefore, the total probability for a given class of 
errors is given by 

p(e) = J2p(se). (5) 

ses 

One can choose a set of undetectable errors Di and use 
them to label the error classes compatible with any given 
syndrome. Namely, if E is compatible with the syndrome 
then the possible error classes are E itself and the classes 
~D~E. 

The error-correction process starts with the measure- 
ment of the check operators Si . Measuring each Si yields 
an eigenvalue Si = ±1. Only certain errors are compati- 
ble with these eigenvalues. In particular, E is compatible 
with the error syndrome if ESi = SiSiE. Ideally, given a 
syndrome s = {s^ one can compute the relative proba- 
bilities P(E\s) of the different error classes E compatible 
with s. If E s is the class that maximizes this probability, 
the best guess is that this is the error that occurred and 
thus should be corrected. The net effect of such an ideal 
error correction is 

Pi > P2 = P0 P + Yl Pi D iP D i ' ( 6 ) 




FIG. 3: (Color online) For the hexagonal arrangement, there 
is a stabilizer operator Z 6 for each of the hexagon plaquettes 
(top left). In the mapping, these stabilizer operators translate 
to classical Ising spins, which are placed on the dual lattice 
(regular triangular lattice, top right). The square-octagonal 
setup (bottom left) has wider computing capabilities because 
it allows for a larger class of quantum gates to be imple- 
mented. There are stabilizers Z 4 [Z 8 ] on the rectangles [oc- 
tagons] . The corresponding dual lattice in the mapping is the 
union jack lattice (bottom right). 



where the success probability po and the probability for 
an effective error Di are 

po := P (Es) Pi := Y, P(^El) ■ (7) 

S S 

Note that in Eq. (j7]) the sum is over possible syndromes. 
Furthermore, 

where D is the number of error classes per error syn- 
drome. In practice this ideal error correction might be 
too costly from a computational perspective. Therefore, 
approximations are needed. 



B. Toric codes and color codes 

Topological codes have two interesting features: First, 
they can be defined for different system sizes in such a 
way that check operators remain local — involving only a 
few neighboring qubits — and at the same time nontrivial 
undetectable errors are global and thus involve a num- 
ber of qubits that depend on the system size. Second, 
they exhibit an error threshold. For error rates below the 
threshold, the success probability [Eq. flTJ)] approaches 1 



for increasing system size, whereas I — po decreases ex- 
ponentially. 

In toric codes [9| physical qubits are placed on the 
edges of a square lattice. Notice that for each edge in 
the direct lattice there is an edge in the dual lattice. 
Check operators Sf are attached to faces /, either in 
the direct or the dual lattices. Toric codes can thus be 
defined in two similar, but distinct ways: In the original 
definition by Kitaev, if / is a face in the direct [dual] 
lattice composed by the edges r, s, t and u, then the 
corresponding check operator is Sf := erf cr^ erf (g) of. 
[Sf := afjg) (T z s af <g) a^\. The second definition is due 
to Wen J27[ and it does not distinguish between dual and 
direct faces. If / has a top edge r, a bottom edge s and 
side edges t,u, then we take Sf := erf (g)crf(g)crf ®cr*. Both 
definitions are equivalent up to a rotation of half of the 
qubits. However, for the depolarizing channel Kitaev's 
definition is related to the alternating eight- vertex model 
and Wen's definition to the standard eight- vertex model. 

In color codes [10[ physical qubits are placed on the 
vertices of a trivalent lattice with three-colorable faces, 
such as, for example, the honeycomb lattice. There are 
two check operators Sf and Sf attached to each face / 
taking the form Sf := i erf and Sf := ®i a fS respec- 
tively, with i running over the qubits on the vertices of 

Because the computing capabilities of color codes de- 
pend on the underlying lattice where the qubits are 
placed, we study two different scenarios: the honey- 
comb lattice for its simplicity, and a lattice of octagons 
and squares that allows for the implementation of ad- 
ditional types of quantum gates. In the mapping onto a 
statistical- mechanical model to compute the error thresh- 
old these two arrangements correspond to the triangular 
and union jack lattices, respectively (see Fig. [3|). 



III. RANDOM EIGHT- VERTEX ISING MODELS 

To determine the error threshold, we show that topo- 
logical codes under the depolarizing channel connect to 
certain random classical spin models. 

For the toric code, the error-correction process maps 
onto a statistical-mechanical interacting eight-vertex 
model [l8l-[2l|. Remarkably, this class of models exhibit 
critical exponents that depend on the coupling constants, 
Eq. ([22]), thus challenging the very notion of universal- 
ity. Eight- vertex models were originally formulated in the 
'electric picture' where the degrees of freedom are elec- 
tric dipoles placed at the bonds surrounding each vertex 
of a square lattice [22[, i.e., the number of independent 
dipole configurations per vertex is eight. In addition, a 
mapping to a 'magnetic picture' was found by Wu [281 ], 
as well as Kadanoff and Wegner (29[: Consider two in- 
dependent Ising systems, each on a square lattice, with 
classical spin variables Si and s' k taking on values ±1, and 
bonds Jij and J' Ml respectively. The lattices are stacked 
as shown in Fig. [2] such that the vertices (spin sites) of 



one lattice are at the center of the plaquettes of the other. 
The Hamiltonian takes the explicit form 

H = - ^(JijSiSj + J' M s f k s f t + J+SiSjs'kSi) . (9) 
+ 

This can be thought of as two interacting Ising models by 
means of a four-spin interaction (denoted by the symbol 
+) between original and dual lattices. 

In fact, two types of eight- vertex models can be related 
to error correction in toric codes: the standard eight- 
vertex model where = J [J^ = J') if a bond is a 
horizontal [vertical] link, see Figs.[2ja) and[2fb); and the 
alternating eight- vertex model where = J [J^ = J') if 
a bond belongs to the direct [dual] lattice. In both cases, 
we set the four-spin interaction to J + = J4, as depicted 
in Fig. [2tc). Thus, both types of eight- vertex models 
share the same lattice structure but differ in the pattern 
of coupling constants. This difference has fundamental 
consequences in the exact solvability of the model: while 
the standard eight- vertex model is exactly solved for ar- 
bitrary couplings, the alternating eight- vertex model is 
generally not exactly solvable. In fact, the latter corre- 
sponds to the Ashkin-Teller model [30] . Notice that when 
J4 = 0, the eight- vertex model reduces to two decoupled 
Ising models, while for J4 7^ the model has two critical 
temperatures. 

The error threshold for correction in quantum codes 
corresponds to the critical line separating ordered from 
disordered phases. The former represents a situation 
where quantum error correction can be performed with 
arbitrary precision. Determining the location of this crit- 
ical line in eight- vertex models is facilitated by the exis- 
tence of a self-duality symmetry in the partition function: 
a duality transformation relating a high-temperature 
eight- vertex model to a low-temperature one on the same 
lattice. Self-duality implies that the coupling constants 
for 2-spin interactions are isotropic, i.e., J = J' . Al- 
together, an isotropic self-dual eight- vertex model has a 
critical line given by [22| : 

J' = J , e~ 2f3j4 = sinh(2/3 J), (10) 

with the restriction that J4 < J. The point in the plane 
( J 4 = J c , J = J c ) at which the self-dual line ceases to be 
critical is given by 

/3J c = ilog(3) « 0.2746.... (11) 

This is already a remarkable and encouraging result be- 
cause it yields a critical point which is approximately 
60% larger than in the standard square-lattice two- 
dimensional Ising model. Note that the error thresh- 
old for bit-flip or phase-flip errors in the Kitaev model 
is computed via a mapping to the aforementioned two- 
dimensional Ising model. In that case, the critical point 
can be computed from the relationship sinh(2/3J c ) = 1, 
i.e., f3J c = 0.4406. Recall that the critical exponents 
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depend continuously on the value of J4. 

In this work we extend the standard eight- vertex model 
by adding quenched disorder to the couplings between 
the spins. Given that for the eight-vertex model f3J c « 
0.2746 is smaller than for the square-lattice Ising model, 
we can expect to find a larger error threshold for the 
depolarizing channel than for bit-flip or phase-flip errors. 

In addition to depolarizing errors in the toric code, 
where the problem map onto Eq. [9j we also study color 
codes, see Fig. [3l In this case the underlying statistical- 
mechanical model to study the error stability to depo- 
larizing errors is defined on a triangular lattice. There 
are two Ising variables — Si and s- — per site. For conve- 
nience, we introduce an artificial third variable s" = s^. 
The Hamiltonian is then given by: 



H 2 = ~ (Jijk SiSjS k + J[ jk 



jii ff ff ff\ 

J ijk S i S j S k ) 



(hj,k) 



(12) 

Equation (fT2|) is illustrated in Fig. 2] where the top [bot- 
tom] layer corresponds to the si [s^] Ising variables with 
the corresponding three-spin interaction term as shown 
in Fig. [4(a) [1(b)]. The third term in the Hamiltonian 
with six-spin interactions is represented by Fig. 2(c). 

When J-j k = in Eq. (fT2]h we obtain two independent 
triangular 3-body Ising models. Interestingly, this model 
can be mapped onto a eight- vertex model on a Kagome 
lattice [3l|. Therefore, the color code Hamiltonian H2 in 
Eq. ([T2|) can be thought of as an interacting eight- vertex 
model (or coupled eight- vertex models). In this work we 
consider two different lattice geometries, triangular and 
union jack (see Fig. [3]). 



IV. MAPPING 
A. Spin models for depolarizing noise 

The goal is to relate the stability of a topological sta- 
bilizer code to depolarizing noise to the ordered phase of 
a suitably chosen classical spin model. However, here we 
consider the more general qubit channel shown in Eq. (pQ) . 
This adds transparency to the mapping and reveals the 
differences between Kitaev's and Wen's versions of the 
toric code with respect to error correction. When Eq. (pQ) 
is applied to each qubit in a code, the net result is a chan- 
nel of the form presented in Eq. (jU). In particular, the 
probability for each Pauli error is 



p(E)=(i- P ) n n 



Pu 



1 



p 



(13) 



where n is the total number of qubits and E w the number 
of appearances of a w in the tensor product forming E. 

The classical spin Hamiltonian is constructed as fol- 
lows: 



2. Associate with each single-qubit Pauli operator a 
an interaction term J^s^s^ • • • such that the 
spins Si correspond to the check operators S l af- 
fected by (j, i.e., such that Sid = —crSi. 

3. Attach to each coupling a sign r a = ±1 dictated by 
the Pauli error E labeling the Hamiltonian, through 
the conditions crE = r a Ea. 

The resulting Hamiltonian takes the general form 

H E = -J2 J - T - s i s 2---SN a , (14) 

cr 

where the sum is over all Pauli operators a and there are 
only three different couplings J a since we set J a w • — J uj 5 
with w = x,y,z and k the qubit label. 

For the mapping, we require the interaction constants 
to be 

T 1 , PxPyPz 

Jw = _ 7^ lo S^71 w w = x,y,z. (15) 

4 £ K0--P) 

This relates the error probability in Eq. ([T3]) to the Boltz- 
mann factor for the ordered state, {s^ = 1}, given the 
interactions generated by E: 



p(E) ex e" 



-PH E ({ Si =l}) 



(16) 



Note that, to recover Eq. (j3]) just notice that replacing 
the error E — >> E' = S j E is equivalent to considering a 
different spin configuration in the original Hamiltonian: 

H s i E {{si}) = H E {{{\ - 2Sij)8i}) . (17) 

Finally, to complete the mapping, the label E in the 
Hamiltonian must describe quenched randomness. In 
particular, the coupling configuration dictated by E ap- 
pears with probability p(E). Equivalently, this means 
that for every qubit k the probability p{r k1 T k1 T k ) for 
each configuration of coupling signs is given by 

p(+l,+l,+l) = 1-p, (18) 

= Px, 

P(-l,+l,-l) = Py, 
= p z . 

In the case of the depolarizing channel 

Px = Py = Pz = p/3 and J x = J y = J z = J. (19) 

The resulting model has two parameters, p and f3J with 
PJ = 1/T. For low p and T the model orders ferromag- 
netically and along the Nishimori line, 



-4/3J 



g/3 

(i- P y 



(20) 



1. Attach a classical spin S{ to each check operator S l . 



this order is equivalent to the noise being correctable. 
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FIG. 4: (Color online) For topological color codes, qubits 
are arranged on trivalent lattices (hexagonal or square- 
octagonal). These codes are mapped to triangular lattices 
(triangular and union jack, respectively) with plaquette in- 
teractions (a,b) on each layer, as well as six-body interactions 
correlating the two layers (c). 

Indeed, comparing error-class probabilities amounts to 
computing free-energy differences 

P ( D i E ) _ Z DiE(P) _. —(3Ai(j3,E) / 91 \ 

P(E) ~ Z E {fi) ' [2l) 

In topological codes we expect the existence of an 
error threshold p c — or several for different error types, 
but we do not need such generality. When p < p c in 
the limit of large systems the success probability is ex- 
pected to approach unity, i.e., po — > 1 and thus due 
to Eq. (|8]) along the Nishimori line the free-energy dif- 
ference is asymptotically infinite, because for any real 
t, P(Ai(P,E) > i) — > 1. Similarly, when p > p c , 
the success probability is expected to become minimal 
(po — > V^) an d thus the free-energy difference con- 
verges in probability zero, so that for any t > we have 
P(\Ai(P,E)\ <t) — y 1. This shows that the-free energy 
differences A$ are order parameters and p c is the criti- 
cal value of p along the Nishimori line. In the models of 
interest here, these are domain- wall free energies. 



B. Models for toric and color codes 

Let us now study what the above mapping, Eq. ([H]h 
produces when applied to toric codes and to topological 
color codes. 

For the toric code, the single-qubit operators <j x and o z 
produce 2-body interactions, because each bit flip [phase 
flip] affects the stabilizer operators on two on two neigh- 
boring dual [direct] faces. The <j y operators, which com- 
bine correlated spin-flip and phase-flip errors, introduce 
four-body interactions, see Fig. [2j The result is an alter- 



nating eight-vertex model with coupling signs that are 
parametrized by a Pauli error E, namely, 

H E = - ^(J x T% S iSj + J z T Z + s' k s[ + J y T\ SiSjS'^) . 

(22) 

The classical spin variables Si [s' k ] live on the top [bot- 
tom] layer of two stacked two-dimensional Ising square 
lattices with interaction constants J x [J z ] (see Fig. [2j). 
The two layers are shifted by half a lattice spacing and 
the third term in Eq. (|22|) describes the combined four 
spin interaction at each of the crossings . Note that 
in Eq. ([22]) there is one qubit per cross +. For Wen's 
toric code one recovers the standard eight- vertex model. 
In either case, there is a global Z2 x Z2 symmetry because 
one can flip all spins in each lattice separately without 
affecting the total energy. 

In the case of color codes there is one spin per face. 
The a x and a z single-qubit operators produce 3-body 
interactions in Eq. ([H|h whereas <j y operators produce 
6-body interactions. The Hamiltonian is then given by 
Eq. ([T2]) but with coupling signs parametrized by a Pauli 
error E, namely, 

He = ~ E Jr»T? jk 8?S?S%, (23) 

(i,j,k) w=x,y,z 

with sfs y s z = 1. Therefore, we obtain two stacked tri- 
angular lattices having three and six-body interactions 
(see Fig. with the six-body interactions introducing 
correlations between both layers. In this case the global 
symmetry is Z2 x Z2 x Z2 x Z2. Indeed, the sites can be 
colored with three colors in such a way that each triangles 
has a site of each color. Thus one can flip all spins for 
two given colors in either of the two lattices separately 
without affecting the total energy. 

For p = in Eqs. (|22|) and ([23]) self-duality predicts a 
critical temperature of T c = 1//3J C ~ 3.641, a value that 
we confirm numerically in our Monte Carlo simulations. 



V. MONTE CARLO SIMULATIONS 

We investigate the classical statistical spin models ac- 
quired in the mapping, Eq. (f22|) and Eq. ([23]) , via large- 
scale classical Monte Carlo simulations using the parallel 
tempering Monte Carlo technique [321 ]. 

In the parallel tempering Monte Carlo method, sev- 
eral identical copies of the system at different tempera- 
tures are simulated. In addition to local simple Monte 
Carlo (Metropolis) spin updates [33[ , one performs global 
moves in which the temperatures of two neighboring 
copies are exchanged. It is important to select the po- 
sition of the individual temperatures carefully such that 
the acceptance probabilities for the global moves are large 
enough J34[ and each copy performs a random walk in 
temperature space. This, in turn, allows each copy to 
efficiently sample the rough energy landscape, therefore 
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speeding up the simulation enormously. 

Detecting the transition temperature T c (p) for differ- 
ent fixed amounts of disorder allows us to pinpoint the 
phase boundary in the p — T phase diagram. The error 
threshold p c is then given by the intersection of the phase 
boundary with the Nishimori line. 



A. Observables and Simulation Details 

For the toric code, it is expedient to partition the lat- 
tice into two sublattices such that the only interconnec- 
tion is given by the four-body-interactions of the Hamil- 
tonian in Eq. ([22]) . The ground state of the pure system 
is realized when the spins of each sublattice are aligned 
(but the alignment may be different as the sign would 
cancel out in both the two and four-spin terms). In this 
case the sublattice magnetization is a good order param- 
eter, 



m 



1 



(24) 



where V denotes one of the sublattices. Similarly, we note 
that each layer of the triangular lattice for color codes is 
tripartite with spins aligned in each sublattice for all re- 
alizations of the pure system's ground state. Hence, we 
can define an order parameter analogous to Eq. ([M]) for 
a suitable sublattice V' . Note that particular caution is 
required when implementing the periodic boundary con- 
ditions to ensure that these distinct sublattices are well 
defined. We can now use the magnetization defined in 
Eq. ([24]) to construct the wave-vector-dependent mag- 
netic susceptibility, 





(25) 



where (• • • )t denotes a thermal average and is the 
spatial location of the spin Si. From Eq. ([25]) we con- 
struct the two-point finite-size correlation function [35| , 



1 



_ / [Xm(0)]a 

2sin(/c min /2) V [Xm(k min )]av 



1, 



(26) 



where [•••] a v denotes an average over disorder and 
kmin = (2ir/L, 0,0) is the smallest nonzero wave vector. 
Near the transition, £c is expected to scale as 



Z L /L~X\l}t v {T-T c . 



(27) 



where X is a dimensionless scaling function. Because 
at the transition temperature, T = T c , the argument 
of Eq. ([27]) becomes zero (and hence independent of L), 
we expect lines of different system sizes to cross at this 
point. If however the lines do not meet, we know that no 
transition occurs in the studied temperature range. 



TABLE I: Simulation parameters: L is the linear system size, 
N sa is the number of disorder samples, t eq = 2 b is the number 
of equilibration sweeps (system size times number of single- 
spin Monte Carlo updates), T m in [T max ] is the lowest [highest] 
temperature, and Nt the number of temperatures used. For 
the toric code, we use L = {12,16,18,24,32,36}, while for 
color codes L = {12, 15, 18, 24, 30, 36} following the coloring 
constraints that the system size must be a multiple of 3. 
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In all simulations, equilibration is tested using a log- 
arithmic binning of the data: This is done by verify- 
ing that all observables averaged over logarithmically- 
increasing amounts of Monte Carlo time are, on average, 
time independent. Once the data for all observables agree 
for three logarithmic bins we deem the Monte Carlo sim- 
ulation for that system size to be in thermal equilibrium. 
The simulation parameters can be found in table |T] 



B. Sample results 

For the pure system (p = 0) there is a sharp transi- 
tion visible directly in the sublattice magnetization. The 
transition temperature T C)Pure w 3.64 coincides with the 
value found analytically. For larger amounts of disorder, 
a transition can still be located precisely by means of 
the crossings in the two-point finite-size correlation func- 
tion [Eq. ([26]) ] for different system sizes. Sample data for 
a disorder strength of p = 0.170 (i.e., this would mean 
that on average 17% of the physical qubits are "broken") 
are shown in Fig. [5j indicating a transition temperature 
of T c (p) = 2.14(2). The error bars are calculated using 
a bootstrap analysis of 500 samples. There are small 
finite-size effects which are addressed by analyzing the 
intersection T*(L, 2L) of pairs of system sizes. We esti- 
mate the limit value for L — >> oo by means of a linear fit 
in a 1 / L-plot - this is our estimate for the best value in 
the physically-relevant thermodynamic limit. For disor- 
der rates approaching the error threshold, corrections to 
scaling increase and a careful finite-size scaling analysis 
has to be performed to determine T c [36|. At p = 0.189, 
the lines only touch marginally such that both the sce- 
nario of a crossing as well as no transition are compatible 
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1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 T 

FIG. 5: (Color online) Crossing of the correlation function 
£,l/L for the toric code with a disorder rate of p = 0.170. The 
data exhibit a clear crossing at a transition temperature of 
T c (p) & 2.14(2). Corrections to scaling are still minimal at 
this disorder rate, but increase closer to the error threshold. 
Inset: closeup of the area where the crossing occurs. The con- 
servative estimate for the transition temperature is indicated 
by the vertical shade. 



within error bars. This gives rise to the large error bars 
in the phase diagram (Fig. [I]). For error rates p > p c , the 
lines do not meet, indicating that there is no transition 
in the temperature range studied. 

The crossing of the critical line T c (p) with the Nishi- 
mori line [Eq. ([20]) ] determines the error threshold to 
depolarization. Our (conservative) estimate is p c = 
0.189(3). Our results are summarized in Fig. [TJ which 
shows the estimated phase diagram. 



VI. DUALITY METHOD 

An alternative approach to estimate the critical value 
p c is to use the duality method (37J , originally developed 
within the context of spin glasses. 

The critical point of a statistical model expressed only 
by local interactions between degrees of freedom can be 
analyzed using the duality method under the assumption 
of a unique transition temperature. The partition func- 
tion Z[A] is then given by the local Boltzmann factor 
= exp(/3Jcos7T(/>), where <p G {0,1} is the difference 
between adjacent spins such as cos(7r</>) = ±1. We de- 
fine the principal Boltzmann factor Aq as the case where 
all spins are parallel. The partition function has to be 
invariant under a Fourier transform, i.e., Z[A] = Z[A*], 
where A* is a dual principal Boltzmann factor (via a 
Fourier transformation). In that case the critical point is 
determined via the equality Aq = Aq. This implies that 
all the components of the local Boltzmann factors are 
equal for several self-dual models such as the standard 
Ising model. Although this self-duality does not work a 
priori for systems with quenched disorder in the general 
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FIG. 6: (Color online) Clusters used to estimate the error 
threshold for the depolarizing channel. The blue lines and 
triangles denote quenched random variables r z , and the red 
lines and triangles correspond to r x . The central site is the 
spin variables summed over. The outer sites represent the 
spin variables fixed in the up direction. 



case, the method can be applied in a special subspace 
called the Nishimori line [37|. The results can be im- 
proved by considering extended local Boltzmann factors 
over a restricted range of interactions [37], HH (see Fig. [6] 
which illustrates the used clusters). Because the result- 
ing statistical-mechanical Hamiltonians for both the toric 
code and topological color codes are self dual, we can ap- 
ply this efficient technique to obtain estimates (up to 
systematic deviations that depend on the clusters used) 
of the error threshold. 



A. Zeroth-order approximation 

The effects of the depolarizing channel on topological 
codes can be expressed by a spin-glass model with the 
partition function [39| 



^] = E£lRi?i 

<f>i 4>i W> 



(28) 



where 



t(rV*) 



- exjp{f3Jr x cos 7r</>+ 
f3Jr z cos 7T(j) + f3Jr x r z cos ircj) cos 7T(/>} . 



(29) 



rfj G {±1} and r?- G {±1} are quenched random vari- 
ables chosen from the distribution 



P(T*,T*0oceW+^ 



(30) 



This model has a gauge symmetry in the subspace J = J p 
which corresponds to the Nishimori line. 

To determine the multicritical point we replicate the 
partition function to take into account the quenched ran- 



domness of the variables rf- and r^-, i.e., 



Z„[/l] 



EE IK 

_ \ <t>i 4>i (ij) 



(t5»t&) 



j,4>i~ 



(31) 



where the brackets denote a configurational average. The 
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local Boltzmann factor is then given by 



u4 



(32) 



where h distinguishes the specific configuration (0f , (f)f ). 
The n-binary Fourier transformation gives the dual 
Boltzmann factor A* k . It follows [13, [HI that A n ^ = 
A* determines the critical point along the Nishimori 
line. Taking the leading term in n, we obtain the er- 
ror threshold for the depolarizing channel of the toric 
code as p c = 0.189... under the conditions J = J p 
and 3exp(— 4/3 J) = p/(l — p) for the Nishimori line. 
Because the local Boltzmann factors for the topological 
color codes on both the hexagonal and square-octagonal 
lattice are the same, we obtain the same estimate for the 
error threshold. 



B. First-order approximation using finite clusters 

To reduce systematic errors we consider finite-size clus- 
ters with four bonds on each square lattice for the toric 
code, six triangles taken from each triangular lattice for 
the color codes on the hexagonal lattice and four trian- 
gles from each union jack lattice for color codes on the 
square-octagonal lattice (see Fig. [6]). We compute the 
principal Boltzmann factors on the clusters, i.e., 



magnetic phase in a random eight-vertex model. Simi- 
larly, color codes turn out to be related to a class of 'in- 
teracting' eight-vertex models. We analyze the models 
resulting from the mapping via both large-scale parallel 
tempering Monte Carlo simulations [16|, [36| and the du- 
ality method [37], By determining T c (p) for different 
error probabilities p, we are able to determine the phase 
boundary in the p — T plane (Fig. [I]). Both approaches 
confirm the existence of a stable ordered phase and by 
locating the intersection of the phase boundary with the 
Nishimori line, we compute in a nonperturbative way, the 
disturbing effects of noise on topological codes. The ex- 
ternal noise considered in this work is more realistic than 
in previous studies because it applies to both bit-flip er- 
rors, phase-flip errors and more importantly, a nontrivial 
combination thereof. 

The error threshold to depolarization errors for differ- 
ent classes of topological codes studied is approximately 
19%, which is larger than the threshold for noncorrelated 
errors. This is very encouraging and shows that topologi- 
cal codes are more resilient to depolarization effects than 
previously thought. The profound relationship between 
complex statistical- mechanical models, such as the eight- 
vertex model, and topological quantum error correction 
promises to deliver a plethora of new paradigms to be 
studied in both fields in coming years. 



4(1) 



Z-^i 11 00-00 
^00,00 W) 

*(1) 



(33) 



as well as its dual A n via a n-binary Fourier transfor- 
mation. As before, the critical point along the Nishimori 

Taking the leading 



line is determined via A 



(i) 

n,0 



order in n we obtain for the error thresholds 

p c = 0.1888... (toric code), (34) 
p c = 0.1914 . . . (color code — hexagonal), (35) 
p c = 0.1878 . . . (color code — sq. octogonal). (36) 

There are small variations in the estimates, however, the 
estimates are all of the order of approximately 19% and 
in agreement with our results from Monte Carlo simula- 
tions. 



VII. RESULTS AND CONCLUSION 

We have shown that the stability under depolarizing 
noise of toric codes can be related to the existence of a 
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